home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / atanh.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  4.9 KB  |  247 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    atanh   double precision hyperbolic arc tangent
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    atanh
  29.  *    machine independent routines
  30.  *    math libraries
  31.  *
  32.  *  DESCRIPTION
  33.  *
  34.  *    Returns double precision hyperbolic arc tangent of double precision
  35.  *    floating point number.
  36.  *
  37.  *  USAGE
  38.  *
  39.  *    double atanh (x)
  40.  *    double x;
  41.  *
  42.  *  RESTRICTIONS
  43.  *
  44.  *    The range of the atanh function is -1.0 to +1.0 exclusive.
  45.  *    Certain pathological cases near 1.0 may cause overflow
  46.  *    in evaluation of 1+x / 1-x, depending upon machine exponent
  47.  *    range and mantissa precision.
  48.  *
  49.  *    For precision information refer to documentation of the
  50.  *    other floating point library routines called.
  51.  *    
  52.  *  PROGRAMMER
  53.  *
  54.  *    Fred Fish
  55.  *    68881 support added by Michael Ritzert
  56.  *
  57.  *  INTERNALS
  58.  *
  59.  *    Computes atanh(x) from:
  60.  *
  61.  *        1.    If x <= -1.0 or x >= 1.0
  62.  *            then report argument out of
  63.  *            range and return 0.0
  64.  *
  65.  *        2.    atanh(x) = 0.5 * log((1+x)/(1-x))
  66.  *
  67.  */
  68.  
  69. #include <stdio.h>
  70. #include <math.h>
  71. #include "pml.h"
  72.  
  73. #if !defined (__M68881__) && !defined (sfp004) /* mjr++        */
  74.  
  75. static char funcname[] = "atanh";
  76.  
  77.     struct exception xcpt;
  78.  
  79. double atanh (x)
  80. double x;
  81. {
  82.     if (x <= -1.0 || x >= 1.0) {
  83.     xcpt.type = DOMAIN;
  84.     xcpt.name = funcname;
  85.     xcpt.arg1 = x;
  86.     if (!matherr (&xcpt)) {
  87.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  88.         errno = ERANGE;
  89.         xcpt.retval = HUGE_VAL;
  90.     }
  91.     } else {
  92.     xcpt.retval = 0.5 * log ((1+x) / (1-x));
  93.     }
  94.     return (xcpt.retval);
  95. }
  96. #endif
  97. #ifdef    sfp004
  98.  
  99. __asm("
  100. comm =     -6
  101. resp =    -16
  102. zahl =      0
  103. ");    /* end asm    */
  104.  
  105. #endif    sfp004
  106. #if defined (__M68881__) || defined (sfp004)
  107.     __asm(".text; .even");
  108.  
  109. # ifdef    ERROR_CHECK
  110.  
  111.     __asm("
  112.  
  113. _Overflow:
  114.     .ascii \"OVERFLOW\\0\"
  115. _Domain:
  116.     .ascii \"DOMAIN\\0\"
  117. _Error_String:
  118.     .ascii \"atanh: %s error\\n\\0\"
  119. .even
  120. | pml compatible atanhgent
  121. | m.ritzert 7.12.1991
  122. | ritzert@dfg.dbp.de
  123. |
  124. |    /* NAN  = {7fffffff,ffffffff}        */
  125. |    /* +Inf = {7ff00000,00000000}        */
  126. |    /* -Inf = {fff00000,00000000}        */
  127. |    /* MAX_D= {7fee42d1,30773b76}        */
  128. |    /* MIN_D= {ffee42d1,30773b76}        */
  129.  
  130. .even
  131. double_max:
  132.     .long    0x7fee42d1
  133.     .long    0x30273b76
  134. double_min:
  135.     .long    0xffee42d1
  136.     .long    0x30273b76
  137. NaN:
  138.     .long    0x7fffffff
  139.     .long    0xffffffff
  140. p_Inf:
  141.     .long    0x7ff00000
  142.     .long    0x00000000
  143. m_Inf:
  144.     .long    0xfff00000
  145.     .long    0x00000000
  146. ");
  147. # endif    ERROR_CHECK
  148.  
  149. __asm("
  150. .even
  151.     .globl _atanh
  152. _atanh:
  153.     ");    /* end asm    */
  154.  
  155. #endif    /* __M68881__ || sfp004    */
  156. #ifdef    __M68881__
  157.  
  158.     __asm("
  159.     fatanhd    a7@(4), fp0    | atanh
  160.     fmoved    fp0,a7@-    | push result
  161.     moveml    a7@+,d0-d1    | return_value
  162.     ");    /* end asm    */
  163.  
  164. #endif    __M68881__
  165. #ifdef    sfp004
  166.     __asm("
  167.     lea    0xfffa50,a0
  168.     movew    #0x540d,a0@(comm)    | specify function
  169.     cmpiw    #0x8900,a0@(resp)    | check
  170.     movel    a7@(4),a0@        | load arg_hi
  171.     movel    a7@(8),a0@        | load arg_low
  172.     movew    #0x7400,a0@(comm)    | result to d0
  173.     .long    0x0c688900, 0xfff067f8    | wait
  174.     movel    a0@,d0
  175.     movel    a0@,d1
  176.     ");    /* end asm    */
  177.  
  178. #endif    sfp004
  179. #if defined (__M68881__) || defined (sfp004)
  180. # ifdef    ERROR_CHECK
  181.     __asm("
  182.     lea    double_max,a0    |
  183.     swap    d0        | exponent into lower word
  184.     cmpw    a0@(16),d0    | == NaN ?
  185.     beq    error_nan    |
  186.     cmpw    a0@(24),d0    | == + Infinity ?
  187.     beq    error_plus    |
  188.     cmpw    a0@(32),d0    | == - Infinity ?
  189.     beq    error_minus    |
  190.     swap    d0        | result ok,
  191.     rts            | restore d0
  192. ");
  193. #ifndef    __MSHORT__
  194. __asm("
  195. error_minus:
  196.     swap    d0
  197.     moveml    d0-d1,a7@-
  198.     movel    #63,_errno    | errno = ERANGE
  199.     pea    _Overflow    | for printf
  200.     bra    error_exit    |
  201. error_plus:
  202.     swap    d0
  203.     moveml    d0-d1,a7@-
  204.     movel    #63,_errno    | NAN => errno = EDOM
  205.     pea    _Overflow    | for printf
  206.     bra    error_exit    |
  207. error_nan:
  208.     moveml    a0@(24),d0-d1    | result = +inf
  209.     moveml    d0-d1,a7@-
  210.     movel    #62,_errno    | NAN => errno = EDOM
  211.     pea    _Domain        | for printf
  212. ");
  213. #else    __MSHORT__
  214. __asm("
  215. error_minus:
  216.     swap    d0
  217.     moveml    d0-d1,a7@-
  218.     movew    #63,_errno    | errno = ERANGE
  219.     pea    _Overflow    | for printf
  220.     bra    error_exit    |
  221. error_plus:
  222.     swap    d0
  223.     moveml    d0-d1,a7@-
  224.     movew    #63,_errno    | NAN => errno = EDOM
  225.     pea    _Overflow    | for printf
  226.     bra    error_exit    |
  227. error_nan:
  228.     moveml    a0@(24),d0-d1    | result = +inf
  229.     moveml    d0-d1,a7@-
  230.     movew    #62,_errno    | NAN => errno = EDOM
  231.     pea    _Domain        | for printf
  232. ");
  233. #endif    __MSHORT__
  234. __asm("
  235. error_exit:
  236.     pea    _Error_String    |
  237.     pea    __iob+52    |
  238.     jbsr    _fprintf    |
  239.     addl    #12,a7        |
  240.     moveml    a7@+,d0-d1
  241.     ");
  242. # endif    ERROR_CHECK
  243.  
  244. __asm("    rts");
  245.  
  246. #endif /* __M68881__ || sfp004    */
  247.